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I.  INTRODUCTION 


The  challenge  to  solve  flow  problems  involving  complicated  geometries  has 
led  to  the  exploration  of  zonal  grid  techniques  using  finite-difference 
methods.  In  the  zonal  grid  approach,  the  global  mesh  for  a  computational  do¬ 
main  is  generated  using  different  zones  separated  by  common  boundaries  called 
zonal  boundaries.  On  the  other  hand,  large  memory  and  substantial  CPU  time 
are  required  for  many  realistic  problems.  On  a  multiprocessor  one  divides  up 
a  program  so  that  many  sections  of  the  code  can  be  worked  on  simultaneously. 
Since  further  gain  in  computational  speed  of  vector  processors  can  only  be 
achieved  at  high  cost  as  upper  limits  on  signal  speed,  packaging  densities  and 
heat  dissipation  are  being  approached,  multiprocessors  evolved  as  higher  level 
architectures.  The  program  can  be  run  much  more  rapidly  if  the  multiprocessor 
can  transmit  data  between  processors,  synchronize  processors  and  make  global 
decisions  efficiently.  Multiprocessor  architectures  can  be  generally  divided 

into  two  categories:  (1)  global  or  shared  memory;  and  (2)  local  or  distribu¬ 
ted  memory.  The  CRAY  X-MP/48  can  be  considered  as  an  example  of  a  shared 

memory  multiprocessor,  while  a  hypercube  can  be  considered  as  an  example  of 
distributed  memory  multiprocessor. 

The  purpose  of  this  report  is  to  document  the  principals  and  operation  of 
a  computer  program  that  has  been  developed  to  solve  the  flow  about  complex 

aerodynamic  shapes  in  a  highly  efficient  manner.  The  computer  code  incorpo¬ 
rates  a  parallel  algorithm  that  can  effectively  utilize  a  distributed  memory 
multiprocessor  architecture  and  a  zonal  grid  approach  that  allows  one  to 
provide  the  capability  to  conveniently  accommodate  complex  boundary  shapes. 

The  code  is  structured  in  such  a  way  that  the  solution  method  is  inde¬ 
pendent  of  zone  coupling  and  parallel  programming  techniques.  This  allows 

execution  of  the  code  on  sequential  or  vector  computers  when  parallel  pro¬ 
gramming  extensions  are  not  invoked.  Also,  incorporation  of  appropriate 
multi-tasking  extensions  of  shared  memory  multiprocessors  allows  the  code  to 
execute  efficiently  on  shared  memory  multiprocessors  such  as  the  CRAY  X-MP/48. 

The  governing  equations  in  general  curvilinear  coordinates  are  discussed 
in  Section  II.  In  Section  III,  a  brief  summary  of  an  algebraic  grid  generator 
is  presented.  The  solution  algorithm,  with  details  of  implementation,  is 
described  in  Section  IV.  Comments  on  the  implementation  of  the  code  and  the 
speed-up  achieved  through  parallel  processing  on  an  INTEL  IPSC  hypercube  and  a 
Cray  X-MP/48  are  given  in  Section  V.  The  computational  code  has  been 
exercised  for  three  flow  field  problems:  constant  and  variable  area  shock 
tube;  tubular  projectile;  and  ramjet  projectile.  Examples  of  computational 
results  are  discussed  and  compared  to  experimental  data  in  Section  VI. 


II.  THEORETICAL  ANALYSIS 

The  compressible,  turbulent  Navier-Stokes  equations  for  axi symmetric  and 
two-dimensional  flow  can  be  expressed  in  the  following  strong  conservation 
form  where  the  dependent  variables  p,  u,  v,  e  are  mass  averaged,  with  e  being 
the  specific  total  energy,  T  the  temperature,  p  and  p  being  mean  density  and 
pressure,  respectively,  and  t  is  time: 
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where  u  is  molecular  viscosity,  e  is  eddy  viscosity  and  M  1  or  0  for  axisym- 
metric  or  two-dimensional  cases,  respectively. 


The  air  was  assumed  to  be  a  perfect  gas  for  both  internal  and  external 
flows,  satisfying  the  equation  state 

p  =  pRT 

where  R  is  the  gas  constant  (1716  ft2/sec2  -  °R  for  air).  For  the  dependence 
of  laminar  viscosity  on  temperature,  Sutherland's  law  was  used: 

v  =  2.270  - 13  --  x  10" 8  .  (2) 

T  +  198.6  ft2 

The  laminar  and  turbulent  Prandtl  numbers,  Pr  and  Prt,  were  assumed  con¬ 
stant  with  values  of  0.72  and  j.9,  respectively.  The  ratio  of  specific  heats, 
y,  was  also  assumed  constant  and  equal  to  1.4.  Cy  and  Cp  are  specific  heat 

capacities  at  constant  volume  and  constant  pressure,  respecti vely. 

(Cy  =  4290  ft2/sec2  -  °R  and  Cp  =  6006  ft2/sec2  -  °R  for  air). 

The  total  energy  per  unit  mass,  e,  is  given  by: 

e  =  CyT  +  (1/2)  (u2  +  v2). 

In  the  5  -  n  computational  plane.  Equation  (la)  is  transformed  to  the 
conservation  law  form  represented  by: 


where 
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The  governing  equations  are  formed  in  such  a  way  that  either  2-D  or  axi- 
symmetric  flow  with  inviscid,  thin  layer  Navi er-Stokes  or  full  Navi er-Stokes 

options  can  be  chosen.  After  some  algebraic  mani pul ati ons  Equation  (3a)  can 
be  transformed  into  the  following  form: 
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The  coefficients  ,  Cj,  (i  =  1,  2. ...4;  j  =  1,2. ...5)  appearing  in  the 
viscous  terms  of  Equation  (4c)  are  defined  as  follows: 
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It  should  be  noted  that  the  viscous  terms  in  Equation  (4a)  have  been 
split  into  terms  which  contain  derivatives  of  flow  variables  in  only  one  di¬ 
rection  (either  £  or  n).  Vectors  and  G2  represent  additional  terms  due  to 

inclusion  of  an  option  for  axisymmetric  2-D  viscous  flow.  When  the  source 
term  (F'  +  H')  and  vectors  and  G2  are  set  to  zero,  one  obtains  a  planar  2-D 

formulation.  It  is  interesting  to  note  that  Equation  (4a)  retains  a  form 
similar  to  that  of  the  Cartesian  counterpart  Equation  (la). 

The  source  term  (F'  +  H')  of  Equation  (4a)  can  be  split  into  inviscid  and 
viscous  terms.  The  viscous  terms  in  turn  can  be  split  into  terms  involving 
derivatives  in  £  and  n  directions. 
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where  AV^j  (i  =  1,2  and  j  =  1,2... .4)  are  defined  as  follows: 
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This  completes  the  derivation  of  the  governing  equations  which  allows  a 
choice  between  2-D  and  axisymmetric  cases  with  inviscid,  thin-layer  or  full 

Navier-Stokes  options. 


III.  GENERATION  OF  COMPUTATIONAL  GRIDS 

As  described  in  Section  II,  the  governing  equations  have  been  expressed 
in  terms  of  general  boundary  conforming  coordinates  in  order  to  facilitate  the 
treatment  of  arbitrary  flow  configurations.  The  technique  employed  in  the 
construction  of  coordinate  systems  is  based  on  algebraic  transformations.  The 
grid  generation  code  reads  in  a  discrete  number  of  user  defined  boundary 
points.  The  code  uses  a  cubic-spline  function  to  define  boundary  shapes  and 
an  analytic  function  to  distribute  the  boundary  points.  The  interior  grid 
points  are  distributed  using  Vinokur's1  clustering  function  based  on  a  hyper¬ 
bolic  tangent.  This  clustering  function  gives  smooth  changes  in  grid  line 
spacing  and  provides  sufficient  control  of  grid  line  spacing  at  the  bound¬ 
aries.  This  grid  generation  procedure  does  not  provide  control  over  grid  line 
orthogonality.  However,  the  procedure  gives  the  desired  degree  of  flexibility 
for  the  problems  under  consideration.  The  non-orthogonal  feature  of  the 
procedure  turns  out  to  be  beneficial  in  that  the  problem  is  not  overly  con¬ 
strained  by  the  orthogonality  condition,  which  would  necessarily  destroy 
boundary  conforming  grids  in  areas  of  sharp  or  highly  concave/convex  corners. 
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Furthermore,  the  grid  distribution  in  the  streamwise  and  radial  directions  are 
independent,  allowing  desired  clustering  at  all  boundaries.  Also,  the 
procedure  permits  the  generation  of  partially  clustered  and  partially  uniform 
grids  in  either  direction.  As  will  be  shown  later,  this  feature  is  useful  for 
some  zonal  grid  overlapping.  Two  examples  which  demonstrate  the  capabilities 
of  this  simple  grid  generation  technique  are  illustrated  in  Figures  1  and  2. 
Grid  clustering  on  the  right  and  upper  boundaries  is  shown  in  Figure  1.  An 
example  where  clustering  is  varied  in  two  regions  of  a  zone  in  the  vertical 
direction  and  a  smooth  variation  in  axial  distribution  achieved  with  different 
left  and  right  increments  is  shown  in  Figure  2. 


IV.  SOLUTION  ALGORITHM 


1.  NUMERICAL  PROCEDURE 


MacCormack 's2  expl icit  and  unsplit  method  is  utilized  for  numerical 
integration  of  the  governing  equations  (4)  in  time  from  an  assumed  initial 
condition  until  a  steady  solution  is  obtained.  The  finite-difference  method 
for  the  one-dimensional  equation: 
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is  given  by  the  following  predictor  -  corrector  steps: 
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where  e"+*  implies  that  the  terms  are  evaluated  using  and  so  forth.  After 

completion  of  the  above  described  two  steps,  first  derivatives  of  the  gov¬ 
erning  equations  are  approximated  by  second-order  accurate  central  differ¬ 
ences.  As  explained  in  Reference  3,  second  derivatives  of  the  viscous  terms 
were  also  effectively  centrally  differenced. 

The  reason  for  using  the  unsplit  method  over  time-split  method  is  to  save 
the  number  of  accessions  of  the  memory.  In  other  words,  for  advancing  one 
time-step,  the  unsplit  method  requires  considerably  less  access  to  the  memory 
than  the  time-split  method.  For  the  explicit  method  the  time-step  size  must 
not  exceed  the  maximum  allowed  by  the  CFL  condition.  An  approximate  linear¬ 
ized  stability  analysis  for  the  inviscid  equations  yields  the  following: 
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where  c  is  the  speed  of  sound.  Since  the  torms  involving  molecular  and  eddy 
viscosity  stabilize  the  solution,  the  time-step  size  computed  using  the  invis- 
cid  analysis  was  found  stable  for  both  inviscid  and  viscous  applications. 
Equation  (8)  was  multiplied  by  a  factor  (denoted  as  CFL)  that  is  slightly  less 
than  one. 

2.  BOUNDARY  AND  INITIAL  CONDITIONS 

The  boundary  conditions  for  the  governing  equations  can  be  categorized 
into  five  major  types  -  freestream,  downstream,  wall,  symmetric  and  no¬ 
reflection.  A  schematic  illustration  of  the  application  of  these  boundary 
conditions  for  a  projectile  flow  field  is  shown  in  Figure  3. 

The  freestream  boundary  conditions  are  held  at  the  appropriate  freestream 
values  for  the  duration  of  the  solution  procedure.  At  the  downstream 
boundary,  the  conventional  zero-gradient  boundary  condition  is  applied.  The 
flow  variables  are  extrapolated  based  on  the  computed  interior  values.  Along 
the  body  surface  n(x,y)  =  0  the  boundary  conditions  are: 


V  =  0  (No  transparency) 

BU 

=  0  (Tangency  or  Slip)  (9) 

=  0  (Adiabatic  Wall) 

while,  for  viscous  flow  the  no-slip  condition  is  U  =  0.  The  pressure  on  the 
body  surface  can  be  obtained  by  applying  normal  pressure  boundary  condition. 

Vp  •  Vn  =  0  (10) 

Using  the  momentum  equations,  the  normal  pressure  boundary  condition  is: 


(Cxnx  +  5yny)p{  ♦  (n|  *  p2)pn  =  p[u{»x)t  *  *(ny)t] 


-  pU(nxu^  +  nyv^) 


(ID 


The  same  relation  has  been  used  in  viscous  flow  with  U  =  0.  the  flow  proper¬ 
ties  along  the  symmetry  line  (n  =  0)  were  specified  by  the  no  gradient  condi¬ 
tions. 
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in  addition  to  v  =  0  being  enforced 
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The  fifth  category  of  boundary  condition  refers  to  the  top  boundary 
(n  =  nmax).  This  boundary  was  established  as  a  no  reflection  boundary.  At 

the  top  boundary,  waves  generated  in  the  interaction  process  in  the  vicinity 
of  the  wall  must  exit  the  computational  domain.  In  order  to  prevent  these 

waves  from  reflecting  from  the  top  boundary  back  into  the  computational  domain 
and  perhaps  contaminating  the  solution,  a  no  reflection  boundary  condition  was 
developed.  The  boundary  conditions  are: 


31 


31  (13) 
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IP  =  0 

31 


where  the  derivative  (3/31)  is  taken  in  the  direction  of  the  outwards  running 
Mach  line.  This  condition  is  appropriate  if  the  flow  is  locally  supersonic 
(i.e. ,  Mach  lines  exist),  and  if  there  is  a  single  family  of  outward  running 
Mach  lines.  A  description  of  this  boundary  condition  is  given  in  more  detail 
in  Appendix  A. 

The  boundary  conditions  are  implemented  in  the  code  using  second-order 
accurate  finite-difference  approximations  in  general. 

For  the  present  calculations,  the  initial  conditions  were  prescribed 
using  either  the  freestream  or  lower  Mach  number  condition  at  all  grid  points 

in  the  computational  domain.  However,  when  available,  previously  computed 
solutions  (obtained  with  different  grids,  etc.)  were  employed  as  initial 
conditions  in  subsequent  calculations. 

3.  NUMERICAL  DAMPING 

Several  types  of  nonlinear  instability  are  encountered  in  the  present 
calculation.  The  remedies  for  these  have  been  discussed  by  MacCormack. 4  One 
of  them  can  be  removed  by  adding  a  fourth-order  pressure  damping  term  to  the 
right  hand  side  of  the  predictor  and  corrector  step.  For  the  S  direction  the 

following  expression  is  added  to  the  flux  j  in  the  predictor  step 


where  a  is  a  damping  constant  with  values  between  0.5  and  5.0.  The  damping 
term  is  of  significance  only  in  the  vicinity  of  pressure  oscillations.  Also, 
this  damping  term  is  more  compact  than  standard  fourth-order  smoothing  terms 

in  that  it  requires  data  at  just  three  grid  points  instead  of  the  usual  five 
points.  It  will  be  shown  later  that,  for  zonal  grids,  this  feature  permits 

the  employment  of  single  grid  cell  overlap. 

4.  ZONAL-OVERLAPPED  GRIDDING 

A  multi-zone  overlapped  gridding  technique  is  used  to  extend  the  time 
dependent  procedure  to  complex  geometries.  The  computational  domain  is 

subdivided  into  several  zones,  each  of  which  require  a  relatively  simple  grid 
generation  technique. 

The  zonal  gridding  approach  has  a  number  of  advantages  over  conventional 
wrap  around  techniques:  (1)  difficulties  in  handling  sharp  corners  can  be 
eliminated,  as  shown  in  Figure  4;  (2)  finer  grids  can  easily  be  used  in 
regions  where  rapid  changes  in  the  flow  variables  occur;  (3)  different  types 

of  equations,  such  as  the  Euler  or  Navi er- Stokes,  can  be  used  in  different 
zones.  (4)  memory  storage  problems  can  be  conveniently  overcome  using 
appropriate  data  management  techniques;  and  (5)  a  saving  in  CPU  time  can  be 
achieved  by  discontinuing  computations  in  the  converged  zones.  The  infor¬ 
mation  exchange  between  zones  must  be  consistent  with  the  governing  equations 
and  lead  to  a  stable  and  efficient  time  marching  or  iterative  scheme.  The 
zonal  boundary  conditions  must  be  transparent  to  shocks  and  regions  of  flow 
separation.  However,  since  zonal  boundary  conditions  introduce  additional 
computational  overhead,  care  must  be  taken  to  insure  that  this  extra  work  is 
significantly  less  than  the  overall  computational  work  required  by  a  solution 
algorithm.  Several  different  techniques  for  zone  coupling  are  discussed  in 
References  5-8. 

A  simple  zone  coupling  technique  is  used  in  the  present  work.  In  this 
technique  zonal  grids  share  one  grid  cell  boundary  with  geometric  continuity 

of  at  least  one  grid  cell  for  overlapped  zones.  The  coupling  of  zones  is 
obtained  by  using  one  grid-cell  overlap.  This  zonal -coupling  is  simple  and 
transparent  to  shock-waves  and  regions  of  flow  separation.  The  transparency 
of  zonal -coupling  is  important  because  the  initial  conditions  are  very  far 
from  the  steady  state  and,  during  the  transient  phase,  shocks  may  travel 
through  the  overlapped  boundaries.  As  shown  in  Figure  5,  two  zones  coincide 
on  a  row  of  overlapped  cells.  The  right-hand  side  boundary  of  zone  A  is 
contained  within  zone  B  and  the  left-hand  side  boundary  of  zone  B  is  contained 
within  Zone  A.  Since  overlapped  cells  are  of  the  same  shape  in  both  zones, 
this  approach  requires  transfer  of  information  from  the  field  of  zone  A  to  the 
boundary  of  zone  B  and  vice-versa. 

Each  multi-zone  solution  was  obtained  by  taking  one  time  step  in  each 
zone  and  then  exchanging  boundary  information  between  zones.  This  zone¬ 
coupling  technique  has  been  found  to  be  very  reliable  and  accurate.  It  is  a 
very  simple  technique  and  well  suited  to  ballistic  projectile  configurations. 
Examples  describing  the  verification  of  the  technique  are  discussed  in  Section 
VI. 
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V.  IMPLEMENTATION  ON  THE  HYPERCUBE  AND  CRAY  X-MP/48 


The  hypercube  is  an  ensemble  of  homogeneous  SISD  (Single  Instruction 
Stream,  Single  Data  Stream)  or  SIMD  (Single  Instruction  Stream,  Multiple  Data 
Stream)  machines  in  a  loosely  coupled  (message  passing),  distributed  (local) 
memory,  MIMD  (Multiple  Instruction  Stream,  Multiple  Data  Stream)  concurrent 
processing  architecture  with  a  hypercube  interconnect  topology.  The  word 
"hypercube"  refers  to  the  specific  way  they  are  interconnected.  A  hypercube 
is  a  generalization  of  the  familiar  three  dimensional  cube,  and  its  topology 
is  completely  specified  by  it's  "dimension."  In  general,  for  a  hypercube  of 

dimension  "d,"  there  are  "2d"  corners,  and  "d"  edges.  Each  corner  of  a  three 
dimensional  cube  is  a  computer,  and  each  edge  is  a  communication  line  between 
the  computers  at  the  two  corners  on  each  end  of  the  edge.  This  is  known  as  a 
three  dimensional  hypercube  (a  "d3").  It  has  eight  "nodes"  (computers),  each 
with  three  nearest  neighbors  and  three  communication  lines  to  those  neighbors. 

The  INTEL  IPSC  hypercube  is  a  "d7"  with  128  nodes  (2^),  and  7  communi¬ 
cation  lines  per  node.  This  hypercube  requires  an  additional  computer  (cube 
manager  or  host)  to  get  programs  and  data  into  the  hypercube  and  results  out. 

A  parallel  algorithm  was  developed  in  such  a  way  that  each  computational 
zone  can  be  assigned  to  a  separate  node  of  a  hypercube.  In  addition,  all  var¬ 
iables  for  each  zone  are  local  so  that  it  can  take  advantage  of  the  local 
memory  architecture.  This  allows  concurrent  computation  in  each  zone  on 
separate  hypercube  nodes.  When  required,  information  can  be  exchanged  between 
overlapped  zones.  This  was  done  using  separate  routines  that  do  zonal  coup¬ 
ling  and/or  synchronization.  Most  of  the  parallel  programming  techniques  for 
local  memory  architectures  such  as  SEND/RECEIVE  are  incorporated  in  these 
routines. 

The  Cray  X-MP/48  is  a  shared  memory  multiprocessor  while  the  hypercube  is 
a  distributed  memory  multiprocessor.  Thus,  different  multi-tasking  techniques 
are  required  when  implementing  an  algorithm  on  shared  or  distributed  memory 
machines.  Again,  mapping  of  zones  onto  processors  is  an  important  issue  for 
efficient  implementation  on  shared  memory  machines.  The  partitioning  across 
processors  for  the  zonal  application  can  be  done  in  varying  degrees  of  granu¬ 
larity.  In  a  coarse  granularity  approach,  one  would  map  a  zone  or  subregion 
of  a  zone  on  each  processor.  In  general,  the  approach  chosen  will  depend  upon 
the  applications.  As  number  of  zones  and/or  size  of  zones  change  from  one 
application  to  another,  the  difficulty  of  uniformly  distributing  work  among 
the  processors  may  require  a  new  mapping  scheme.  This  load  imbalance  would 
result  in  processors  remaining  idle  until  the  processor  with  the  largest  task 
completes  its  portion  of  the  computation.  For  the  eight-zone  ramjet  problem 
each  zone  was  assigned  to  a  processor.  This  was  implemented  using  $  PROCESS 
and  $  ALSO  PROCESS  directives  described  by  Misegades.9  In  going  from  one  to 
four  processors  a  speed-up  of  2.3  was  achieved  on  the  Cray  X-MP.  In  a  medium 
granularity  approach,  all  processors  might  work  on  a  single  zone  in  parallel. 
This  parallelization  consists  of  decomposing  a  zone  to  DO-LOOP  level.  Since 
no  artificial  boundaries  are  introduced  on  shared  memory  machines,  this 
approach  generally  allows  for  better  load  balancing  than  a  coarse  granularity 
approach.  On  the  Cray  this  can  be  implemented  using  microtasking  directives 
such  as  00  GLOBAL.  For  the  eight-zone  ramjet  application,  this  approach 
achieved  speed-up  of  3.55  in  going  from  one  to  four  processors.  This  speed-up 
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could  be  less  than  the  maximum  speed-up  that  can  be  obtained  for  the  applica¬ 
tion  because  several  computations  such  as  boundary  conditions  and  zonal  condi¬ 
tions  were  not  done  in  parallel. 


VI.  COMPUTATIONAL  RESULTS 

In  order  to  evaluate  the  code,  a  series  of  test  problems  were  solved.  It 
was  necessary  to  check  whether  any  perturbation  was  created  in  the  flow  field 
by  the  use  of  the  zonal  boundary  conditions.  Also,  it  was  necessary  to  prove 
that  zonal  boundaries  allow  shocks  to  travel  across  the  interfaces.  This  was 
required  because,  for  many  applications,  the  initial  conditions  are  very  far 
from  the  steady  state  and,  while  the  solution  is  converging,  shocks  travel 
through  the  zonal  interfaces. 

1.  CONSTANT  AND  VARIABLE  AREA  SHOCK  TUBE 

A  simple  constant  area  shock  tube  case  was  considered  for  inviscid  2-D 
and  axisymmetric  flow.  The  purpose  of  this  initial  test  case  was  to  evaluate 
the  ability  of  the  zonal  boundary  technique  to  pass  disturbances  across  zonal 
boundaries  without  creating  pertubations  in  the  flow  field.  The  purpose  was 
not  to  give  detailed  comparison  for  a  specific  pressure  ratio  in  terms  of  over 
pressure,  etc.;  but  to  compare  the  output  of  a  single  zone  case  with  that  of  a 
multi-zone  case  at  every  time-step.  Multi-zone  cases  ranging  from  two  zones 
to  16  zones  were  considered  (Figure  6).  The  results  of  each  multi-zone  test 
case  matched  the  single  zone  case  bit  by  bit.  No  perturbation  was  found  and 
the  shock  traveled  through  the  entire  computational  domain.  Next,  a  variable 
cross-sectional  area  shock  tube  was  considered  as  a  test  case.  The  shock  tube 
had  a  variable  diameter  driver  section,  a  convergent-divergent  nozzle  with  a 
diaphram  at  the  throat  and  a  constant  area  driven  section.  Again,  under  the 
same  flow  conditions,  comparison  between  the  results  for  a  single  zone  and 
multi-zone  cases  (not  shown)  confirmed  the  validity  of  the  zonal  boundary 
technique. 

2.  TURBULAR  PROJECTILE  AT  SUPERSONIC  VELOCITIES 

In  order  to  evaluate  the  axisymmetric  Navier-Stokes  code,  the  flow  in  and 
around  a  tubular  projectile  was  investigated.  The  tubular  configuration  is 
shown  in  Figure  7.  As  shown  in  Figure  8  the  computational  domain  was  sub¬ 
divided  into  five  overlapped  zones.  The  total  grid  of  105  *  30  for  the  inter¬ 
nal  flow  and  95  x  24  for  the  external  flow  was  used  for  obtaining  the  solu¬ 
tions.  Figures  9a  and  9b  show  the  computed  Mach  number  contours  for  inviscid 
and  viscous  solutions,  respectively,  for  the  freestream  Mach  number  of  1.7. 
The  Reynolds  number  is  2.75  x  io6  based  on  chord  length  for  the  laminar 
viscous  case.  A  Schlieren  photograph  for  the  same  case  is  shown  in  Figure  10. 
In  both  figures,  a  strong  shock  (choked  flow)  is  visible  at  the  leading  edge. 
Figures  11  and  12  show  computed  Mach  contours  for  inviscid  and  viscous  solu¬ 
tions,  respectively,  for  M^  =  3.5.  Figure  13  shows  a  Schlieren  photograph  for 

the  freestream  Mach  number  of  3.5  and  Reynolds  number  5.67  x  io6.  For  this 
case,  a  weak  oblique  shock  is  seen  at  the  leading  edge.  The  pressure  distribu¬ 
tion  as  a  function  of  axial  position  along  the  internal  surface  for  both  cases 
are  shown  in  Figures  14  and  15,  respectively.  The  pressure  distribution  is 
about  the  same  for  both  inviscid  and  viscous  flow.  The  difference  in  the 
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pressure  distribution  between  Mach  1.7  and  Mach  3.5  can  be  explained  in  view 
of  the  normal  and  oblique  shock  structures  in  the  forebody.  For  M^  =  3.5,  the 

pressure  rise  at  the  surface,  near  the  exit  of  the  tubular  projectile  is 
interesting.  This  complex  behavior  can  be  explained  with  the  help  of  the  Mach 
contour  plot.  In  the  forebody  region,  an  oblique  shock  generated  by  the  lip 
of  the  tubular  projectile  reflects  off  the  axis.  This  reflecting  shock 
interacts  with  flow  that  is  expanding  over  an  inverted  wedge  corner  which  is 
located  at  the  center  of  the  projectile.  This  shock  boundary  layer  inter¬ 
action  is  believed  to  be  responsible  for  the  pressure  rise  at  the  surface  in 
the  exit  region  of  the  tubular  projectile.  This  problem  was  successfully 
implemented  on  the  INTEL  hypercube  multiprocessor. 

3.  SOLID  FUEL  RAMJET  (SFRJ) 

Finally,  the  flow  in  and  around  a  solid  fuel  ramjet  (SFRJ)  configuration 
was  considered.  The  design  objective  for  the  tubular  ramjet  powered  projec¬ 
tile  is  to  simulate  the  flight  trajectory  of  a  kinetic  energy  (KE)  projectile. 
A  schematic  illustration  of  the  SFRJ  projectile  along  with  the  type  of 
boundary  conditions  imposed  are  shown  in  Figure  16a.  The  internal  configura¬ 
tion  of  the  SFRJ  consists  of  a  diverging  inlet  formed  by  an  axi symmetric  cowl 
followed  by  a  circular  injector  plate  (or  flame  holder),  a  combustor  and  a 
nozzle.  The  combustor  is  located  between  the  injector  plate  and  converging- 
diverging  nozzle.  Solid  fuel  is  located  along  the  circular  wall  of  the 
combustor.  The  external  configuration  consists  of  a  tangent-ogive  forebody 
followed  by  a  straight  circular  afterbody. 

A  series  of  wind  tunnel  tests10  were  initiated  in  order  to  provide  data 
to  help  guide  the  development  of  computational  modeling  techniques  for  the 
SFRJ  configuration  and  to  provide  guidance  for  the  design  of  a  prototype 
projectile.  The  wind  tunnel  test  conditions  were:  freestream  Mach  Number  = 
4.0;  freestream  stagnation  pressure  =  14.36  psia;  and  freestream  stagnation 
temperature  =  480  R.  The  experimental  data  consisted  of  internal  surface 
pressure  measurements  and  flow  visualization. 

The  objectives  of  this  computational  study  were  to:  (1)  simulate  cold 
flow  inside  a  standard  SFRJ  (nozzle  diameter  1.54  inches  and  injector  diameter 
1.7  inches);  (2)  simulate  a  normal  shock  inlet  condition  by  varying  the  nozzle 
throat  diameter;  and  (3)  study  the  internal  flow  using  graphical  flow  visuali¬ 
zation  techniques  in  order  to  provide  a  better  understanding  of  the  experi¬ 
mentally  determined  internal  surface  pressure  distributions.  Modeling  of  the 
combustion  process  is  of  interest  and  will  be  addressed  in  the  future. 

Because  of  the  complex  internal  geometry  which  involves  several  sharp 
corners,  a  zonal  gridding  approach  was  used.  The  computational  domain  was 
subdivided  into  eight  zones  (Figure  16a).  In  each  zone,  an  algebraic  grid  was 
used.  A  single-cell  overlap  (Section  IV)  was  used  to  couple  neighboring 
zones.  Preliminary  runs  with  an  eight-zone  grid  indicated  sonic  flow  at  the 
nozzle  throat  and  under-expanded  flow  at  the  nozzle  exit.  Since  the  flow  in 
the  diverging  section  of  the  nozzle  remained  supersonic,  the  flow  in  the  base 
area  had  very  little  upstream  influence.  To  save  CPU  time,  a  substantial 
number  of  grid  points  were  removed  from  zone  four  and  zone  seven,  and  zone 
eight  was  totally  eliminated.  Computations  were  performed  for  both  Euler  and 


Navi er-Stokes  cases.  Computational  results  are  *  di scussed  below  for  four 
different  nozzle  throat  diameters  for  an  injector  diameter  of  1.7  inches. 

Figure  16b  shows  a  seven-zone  grid  for  inviscid  computations.  The  grid 
is  almost  uniform  in  the  internal  region.  The  outflow  boundaries  of  zones 
four  and  seven  were  truncated.  The  grid  size  for  zones  one  through  seven  was 
15  x  35,  46  x  35,  15  x  21,  26  x  20,  7  x  25,  84  x  30,  and  12  x  23,  respective¬ 
ly.  Inviscid  computations  for  nozzle  throat  diameters  greater  than  1.45 
showed  an  internal  wall  pressure  level  (P/P.J  of  less  than  four  in  both  the 
cowl  and  combustor  regions.  As  shown  in  Figure  17,  for  a  nozzle  throat 
diameter  of  1.42  inches,  the  pressure  level  in  the  cowl  area  is  close  to  the 
wind  tunnel  measurement  (see  Figure  20)  for  a  nozzle  throat  diameter  of  1.54 
inches.  For  the  wind  tunnel  measurement,  the  equivalent  nozzle  throat 
diameter  is  less  than  1.54  inches  because  of  the  boundary  layer  displacement 
thickness.  Although  the  inviscid  computation  captured  the  oblique  shock  and 
the  pressure  level  in  the  inlet  region  quite  well,  it  failed  to  obtain  the 
measured  pressure  rise  in  the  combustor.  A  Mach  number  contour  plot  (Figure 
17)  shows  an  oblique  shock  originating  just  behind  the  lip  of  the  cowl. 

Computations  for  the  nozzle  throat  diameter  of  1.28  inches  showed  an 
oscillatory  shock  behavior.  Figures  18a  and  18b  show  a  strong  oblique  shock 
becoming  normal  at  the  cowl  lip.  However,  this  normal  shock  does  not  remain 
steady  at  the  cowl  lip.  This  normal  shock  starts  moving  out  of  the  cowl  and 
becomes  a  bow  shock  (Figures  18c  and  18d).  Again  the  bow  shock  shows  unsteady 
behavior  and  it  becomes  a  normal  shock  at  the  lip  and  a  strong  oblique  shock. 
The  above  described  behavior  was  found  cyclic.  In  other  words,  the  flow  keeps 

going  through  a  cycle  of  strong  oblique  to  normal  to  detached  bow  shock  and 
reversal  of  the  sequence.  The  flow  inside  the  cowl  and  combustor  is  subsonic 

for  the  normal  and  detached  bow  shock  conditions.  Existence  of  oscillatory 
entry  shock  phenomena  was  observed  in  aerodynamic  range  tests  of  a  20mm  SFRJ 
(Reference  11). 

Figure  19  shows  a  seven  zone  grid  used  for  the  Navi er-Stokes  computa¬ 
tions.  The  grid  size  for  zones  one  through  seven  was  20  x  64,  46  x  64, 
20  x  25,  64  x  24,  10  x  32,  84  x  52,  and  41  x  32,  respectively.  To  resolve 
high  gradients  near  walls  the  grid  lines  were  clustered  near  walls.  A  laminar 
flow  computation  for  a  1.54  inch  nozzle  indicated  a  constant  wall  pressure 
level  of  about  five  in  the  inlet  and  combustor  sections.  A  Mach  number 
contour  plot  showed  supersonic  core  flow  along  the  axis.  Separated  primary 
and  secondary  flows  were  noticed  inside  the  inlet  and  combustor.  A  shock 
reflection  pattern  consisting  of  four  normal  shocks  at  the  axis  was  observed. 
A  detailed  discussion  of  the  separated  inlet  flow  of  the  present  computations 
can  be  found  in  Reference  12.  For  the  laminar  flow  computation,  the  computed 
wall  pressure  in  the  inlet  section  was  in  fair  agreement  with  data.  However, 
the  computed  wall  pressure  for  the  laminar  flow  did  not  match  the  wind  tunnel 
data  in  the  combustor  section.  This  case,  with  a  nozzle  throat  diameter  of 
1.54  inches,  was  found  to  be  highly  sensitive  to  small  changes  in  model  throat 
diameter.  This  case  is  thought  to  be  sensitive  due  to  the  following  three 
conditions:  (1)  separated  flow  inside  the  combustor;  (2)  boundary  layer 
displacement  thickness  at  the  nozzle  throat;  and  (3)  internal  shock  reflection 
pattern.  A  less  than  satisfactory  prediction  of  these  conditions  was  thought 
to  be  a  potential  reason  for  the  discrepancy  between  computed  and  measured 
data  in  the  combustor.  As  shown  in  Figure  20,  this  highly  sensitive  behavior 
is  visible  in  experimental  measurements  for  three  different  wind  tunnel  runs 
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for  nozzle  throat  diameters  of  1.54  inches  and  1.534  inches  (Reference  13). 
In  an  attempt  to  simulate  the  separated  flow,  a  modified  Bal  dwi  n- Lomax  turbu¬ 
lence  model 14  was  incorporated  in  the  code.  The  computed  wall  pressure  with 
this  turbulence  model  showed  a  trend  that  is  similar  to  one  that  was  observed 
experimentally  inside  the  combustor  (Figure  20).  Although,  this  modified 
model  seems  to  help  for  this  particular  case,  critical  evaluation  based  on 
separated  flow  test  cases  will  be  required  before  making  any  specific 
conclusion  about  the  model.  A  plot  of  Mach  number  contours  (Figure  21)  shows 
a  series  of  shock  reflections  which  coalesce  into  a  normal  shock  in  the  middle 
of  the  combustor.  The  flow  between  this  shock  and  the  nozzle  throat  is 
essentially  subsonic.  A  normalized  pressure  contour  plot  for  the  same  case  is 
shown  in  Figure  22.  To  obtain  a  normal  shock  at  the  lip  of  the  inlet,  a 
series  of  computations  were  made  with  gradual  reduction  in  the  nozzle  throat 
diameter.  At  a  nozzle  throat  diameter  of  1.21  inches,  a  normal  shock  at  the 
lip  of  the  inlet  was  observed  (Figure  23).  The  computed  wall  pressure  level 
shown  in  Figure  24,  is  in  fair  agreement  with  wind  tunnel  data  (not  shown)  for 
a  nozzle  throat  diameter  of  1.25  inches  as  shown  in  Figure  24,  see  Reference 
10. 

CONCLUDING  REMARKS 

A  computer  code  has  been  developed  that  incorporates  a  zonal  grid  techni¬ 
que  to  handle  complex  flow  configurations  and  a  parallel  algorithm  to  exploit 
new  multiprocessor  architectures.  The  code  has  been  exercised  on  the  Intel 
IPSC  hypercube  and  CRAY  X-MP/48  computers  to  verify  the  ability  to  perform 
parallel  computations  for  realistic  fluid  dynamic  problems. 

The  results  of  several  test  cases  which  demonstrated  the  ability  of  the 
zonal  technique  to  compute  highly  complex  flow  cases  have  been  discussed. 

It  is  felt  that  the  parallel  programming  techniques  are  working  very  well 
and  that  the  code  has  excellent  potential  for  performing  computations  for 

complex  configurations  such  as  guided  and  non  axisymmetric  projectiles. 
Further  development  of  the  code  is  continuing  with  emphasis  placed  on  exten¬ 
sions  to  three-dimensions  and  multi-tasking  on  the  CRAY  X-MP/48. 
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Figure  1 .  Fully  clustered  grid  for  cowl  region 
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Figure  2.  Partly  clustered  partly  uniform  grid  for  cowl  region 
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Figure  3.  Zonal  boundary  conditions 


|W‘ 


•  SUBDOMAIN  A 


-  -X  -  -x  - 


i  •  i 

•*--*-* 


--X 


X  SUBDOMAIN  B 


Figure  5.  One  grid-cell  overlap 
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Figure  7.  Experimental  supersonic  tubular  projectile  model 
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Figure  8.  A  five-zone  grid  for  the  tubular  projectile 
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Figure  10.  Wind  tunnel  Schlieren  photograph, 
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tic  illustration  of  ramjet  and  boundary  conditions 
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Figure  16b.  Seven-zone  grid  for  the  ramjet 
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Figure  18a.  Mach  number  contours  showing  a  strong  oblique  shock  in  the  inlet  region 
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Figure  18b.  Mach  number  contours  showing  a  normal  shock  in  the  inlet  region 
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INVISCID 

NOZZLE  DIAMETER  =  1.28 


Figure  18c.  Mach  number  contours  showing  an  attached  normal  shock  at  the  cowl  lip 
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Figure  I8d.  Mach  number  contours  showing  a  detached  bow  shock 
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Figure  20.  Interna!  surface  pressure  distribution  for  the  ramjet  with  1.54  inch  nozzle  diameter 


Figure  22.  Normalized  pressure  contours  for  the  ramjet  with  1 .54  inch  nozzle  diameter 
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Figure  24.  Internal  surface  pressure  distribution 
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LIST  OF  SYMBOLS 


B.j  coefficient  in  Navi er-Stokes  equation 

Cj  coefficient  in  Navi er-Stokes  equation 

c  speed  of  sound 

Cp  specific  heat  at  constant  pressure 

Cv  specific  heat  at  constant  volume 

D.j  coefficient  in  Navi  er-Stokes  equation 

E  flux  vector  in  £  direction 

e  total  energy  per  unit  mass 

F  flux  vector  in  n  direction 

Gj,G2  vectors  of  axi symmetric  viscous  terms 

H  source  term  in  Navi er-Stokes  equation 

J  Jacobian  of  coordinate  transformation 

M  Mach  number 

Pr  molecular  Prandtl  number 

Prt  turbulent  Prandtl  number 

P  static  pressure 

Q  vector  of  dependent  variables 

R  gas  constant 

Sj^  vectors  of  viscous  terms 

T  static  temperature 

^l’^2  vectors  of  viscous  terms 

t  time 

u,v  Cartesian  velocity  components 

U,V  contravariant  velocity  components 

x,y  physical  Cartesian  coordinates 

a  pressure  damping  coefficient 
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LIST  OF  SYMBOLS  (Continued) 


u  molecular  dynamic  viscosity 

e  turbulent  eddy  viscosity 

transformed  coordinates 
p  density 

Subscripts 

i ,j  mesh  indices 

®  free  stream  condition 
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APPENDIX  A:  NON-REFLECTING  OUTER  BOUNDARY  CONDITION 


This  appendix  describes  the  derivation  of  the  nonreflecting  outer 
boundary  condition  employed  along  the  top  boundary  (Figure  3).  The  non¬ 
reflection  boundary  condition,  like  many  other  Neumann  type  boundary  condi¬ 
tions,  is  an  approximation.  For  appropriate  supersonic  flow  applications,  the 
use  of  non-reflection  boundary  conditions  allows  the  outer  boundary  to  be 
placed  close  to  the  body.  This  substantially  reduces  the  number  of  grid 
points  required  for  a  given  application  and  reduces  both  computer  storage  and 
CPU  time  requirements.  This  boundary  condition  is  consistent  when  the  follow¬ 
ing  conditions  are  met. 

a.  The  local  Mach  number  is  greater  than  one  at  the  boundary. 

b.  There  is  a  family  of  outward  running  Mach  lines. 

Consider  a  nonorthogonal  coordinate  system  as  shown  in  Figure  Al. 
Through  point  x,y  one  can  draw  the  velocity  vector  and  Mach  wave  inclined  at 

an  angle  y  with  respect  to  the  velocity  vector.  The  velocity  vector  is  at  an 
angle  measured  with  respect  to  Cartesian  coordinate  x.  This  gives 


M 

and 


id  =  tan-1 


which  leads  to  a  total  included  angle 


0  =  y  +  u) 


let 


f  =  f(p,  u,  v,  T)T. 


The  non-reflection  condition  requires  that  the  flow  variables  remain 
constant  along  the  left  running  Mach  lines. 

fi,JL  =  fx,y 


which  gives 


f 


k-1 


* 


JL-1 


_A _ 

'  fk-l,JL-l) 


where 


x 


(yk,JL-l  ~  yi,JL  4  xi,JL 


tan  0  -  x^ 


,JL-1 


(tan  e  -  tan  Y) 
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tan  f) 


y  .  yjjJL  +  tan  9  (x-x1JL) 


A  =  {(x 


Ax  =  xk,JL-l  "  xk-l,  JL-1 

Ay  =  yk,JL-l  "  yk-l,  JL-1 
¥  =  tan"1  (Ay/ Ax) 

*k-i,  jl-i>2  +  O'  ■  »w,  jl-i)2)1/2 

B  =  (Ax2  +  Ay2) l^z 
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USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  publishes. 
Your  comments/answers  to  the  items/questions  below  will  aid  us  in  our  efforts. 

1.  BRL  Report  Number _ BRL -MR-3834 _ Date  of  Report _ MAY  1990 _ 

2.  Date  Report  Received _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or  other  area  of  interest 

for  which  the  report  will  be  used.) _ 


4.  Specifically,  how  is  the  report  being  used?  (Information  source,  design  data,  procedure,  source 
of  ideas,  etc.) _ 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or  dollars 
saved,  operating  costs  avoided,  or  efficiencies  achieved,  etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate 
changes  to  organization,  technical  content,  format,  etc.)  _ _ 


1,1  rr .'3T 


Name 


CURRENT  Organization 

ADDRESS  _ 

Address 


City,  State,  Zip  Code 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  New  or  Correct 
Address  in  Block  6  above  and  the  Old  or  Incorrect  address  below. 


Name 


OLD  Organization 

ADDRESS  _ 

Address 


City,  State,  Zip  Code 


(Remove  this  sheet,  fold  as  indicated,  staple  or  tape  closed,  and  mail.) 


. FOLD  HERE 

Department  of  the  Army 

Director 

U.S.  Army  Ballistic  Research  Laboratory 

ATTN:  SLCBR-DD-T 

Aberdeen  Proving  Ground,  MD  2U)t'  -5066 
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FIRST  CLASS  PERMIT  No  0001,  APG,  MD 
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U.S.  Army  Ballistic  Research  Laboratory 

ATTN:  SLCBR-DD-T 

Aberdeen  Proving  Ground,  MD  21005-9989 
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